Cargar datos en R, o incluso omitir los datos vacíos (NA), filtrado de la serie (mediante paquetería Tsibble).
Gráfica de los datos (visualización)
global_economy %>%filter(Country =="Sweden") %>%autoplot(GDP) +ggtitle("PIB de Suecia") +ylab("$US billions")
Definición del modelo (especificación)
Describir el modelo. Hay distintos tipos, se usará por ejemplo, un modelo lineal de series de tiempo con TLSM:
TSLM(GDP ~trend())
<TSLM model definition>
Entrenamiento del modelo (estimación)
Una vez se define el modelo, lo que sigue es entrenar al modelo. Lo que significa pasarle los datos para que, estdísticamente, encuentre los parámetros que realizan el mejor ajuste posible.
fit <- global_economy %>%model(Modelo_tendencia =TSLM(GDP ~trend()))
Ahora tenemos un modelo lineal ajustado y el objeto resultante es un mable (model table).
fit
# A mable: 263 x 2
# Key: Country [263]
Country Modelo_tendencia
<fct> <model>
1 Afghanistan <TSLM>
2 Albania <TSLM>
3 Algeria <TSLM>
4 American Samoa <TSLM>
5 Andorra <TSLM>
6 Angola <TSLM>
7 Antigua and Barbuda <TSLM>
8 Arab World <TSLM>
9 Argentina <TSLM>
10 Armenia <TSLM>
# … with 253 more rows
Revisar el desempeño del modelo (evaluación)
Ver como se ajusta el modelo a los datos y compararlo con otros modelos.
Producir pronósticos
Una vez tenemos un modelo ajustado, podemos hacer pronósticos. Usamos forecast(), h = periodo a pronosticar
fcst <- fit %>% fabletools::forecast(h =36) # h = "3 years"fcst
En este pronósitco, las predicciones de todos los valores futuros son la media de los datos históricos.
bricks %>%model(MEAN(Bricks))
# A mable: 1 x 1
`MEAN(Bricks)`
<model>
1 <MEAN>
Método ingenuo (Naive method)
Se toma el último valor como el pronóstico para todos los valores futuros. También se les conoce como pronósitcos de caminata aleatoria, ya que son óptimos en una serie que siga una.
bricks %>%model(NAIVE(Bricks),RW(Bricks)) # hace exactamente lo mismo que NAIVE()
# A mable: 1 x 2
`NAIVE(Bricks)` `RW(Bricks)`
<model> <model>
1 <NAIVE> <NAIVE>
Método ingenuo estacional (seasonal Naive)
Se le agrega un componente a Naive para lidiar con datos altamente estacionales.
bricks %>%model(SNAIVE(Bricks))
# A mable: 1 x 1
`SNAIVE(Bricks)`
<model>
1 <SNAIVE>
Warning: There was 1 warning in `filter()`.
ℹ In argument: `time_in(Time, ...)`.
Caused by warning:
! Argument 'roll' is deprecated. Deprecated in version '1.8.4'.
Método de la deriva (drift)
Variación de Naive que permite que aumente o disminuya el pronótico en el tiempo. El aumento del cambio es el cambio promedio en los datos históricos.
bricks %>%model(RW(Bricks ~drift()))
# A mable: 1 x 1
`RW(Bricks ~ drift())`
<model>
1 <RW w/ drift>
Ejemplos
# Set training data from 1992 to 2006train <- aus_production %>%filter_index("1992 Q1"~"2006 Q4")# Fit the modelsbeer_fit <- train %>%model(Mean =MEAN(Beer),`Naïve`=NAIVE(Beer),`Seasonal naïve`=SNAIVE(Beer),Drift =RW(Beer ~drift()) )# Generate forecasts for 14 quartersbeer_fc <- beer_fit %>%forecast(h=14)# Plot forecasts against actual valuesbeer_fc %>%autoplot(filter_index(aus_production, "1992 Q1"~ .), level =NULL) +ggtitle("Forecasts for quarterly beer production") +xlab("Year") +ylab("Megalitres") +guides(colour=guide_legend(title="Forecast")) +geom_vline(xintercept =as.Date("2007-01-01"), color ="firebrick",linetype ="dashed") +annotate("label", x =c(as.Date("2003-01-01"),as.Date("2009-01-01")),y =550, label =c("Train set", "Test set"),color =c("black","blue"))
beer_fc %>%filter(.model =="Seasonal naïve") %>%autoplot(filter_index(aus_production, "1992 Q1"~ .)) +ggtitle("Seasonal naive forecast for quarterly beer production") +xlab("Quarter") +ylab("Megalitres") +geom_vline(xintercept =as.Date("2007-01-01"), color ="firebrick", linetype ="dashed") +annotate("label", x =c(as.Date("2003-01-01"), as.Date("2009-01-01")), y =550, label =c("Train set", "Test set"), color =c("black", "blue"))
Otro ejemplo:
gafa_stock %>%distinct(Symbol)
# A tibble: 4 × 1
Symbol
<chr>
1 AAPL
2 AMZN
3 FB
4 GOOG
# Re-index based on trading daysgoogle_stock <- gafa_stock %>%filter(Symbol =="GOOG") %>%mutate(day =row_number()) %>%update_tsibble(index = day, regular =TRUE)# Filter the year of interestgoogle_2015 <- google_stock %>%filter(year(Date) ==2015)# Fit the modelsgoogle_fit <- google_2015 %>%model(Mean =MEAN(Close),`Naïve`=NAIVE(Close),Drift =NAIVE(Close ~drift()),SNAIVE =SNAIVE(Close) )
Warning: 1 error encountered for SNAIVE
[1] Non-seasonal model specification provided, use RW() or provide a different lag specification.
google_fit
# A mable: 1 x 5
# Key: Symbol [1]
Symbol Mean Naïve Drift SNAIVE
<chr> <model> <model> <model> <model>
1 GOOG <MEAN> <NAIVE> <RW w/ drift> <NULL model>
# Produce forecasts for the 19 trading days in January 2015google_fc <- google_fit %>%forecast(h =19)google_fc
# A fable: 76 x 5 [1]
# Key: Symbol, .model [4]
Symbol .model day Close .mean
<chr> <chr> <dbl> <dist> <dbl>
1 GOOG Mean 505 N(602, 6766) 602.
2 GOOG Mean 506 N(602, 6766) 602.
3 GOOG Mean 507 N(602, 6766) 602.
4 GOOG Mean 508 N(602, 6766) 602.
5 GOOG Mean 509 N(602, 6766) 602.
6 GOOG Mean 510 N(602, 6766) 602.
7 GOOG Mean 511 N(602, 6766) 602.
8 GOOG Mean 512 N(602, 6766) 602.
9 GOOG Mean 513 N(602, 6766) 602.
10 GOOG Mean 514 N(602, 6766) 602.
# … with 66 more rows
# A better way using a tsibble to determine the forecast horizonsgoogle_jan_2016 <- google_stock %>%filter(yearmonth(Date) ==yearmonth("2016 Jan"))google_fc <- google_fit %>%forecast(google_jan_2016)# Plot the forecastsgoogle_fc %>%autoplot(google_2015, level =NULL) +autolayer(google_jan_2016, Close, color='black') +ggtitle("Google stock (daily ending 31 Dec 2015)") +xlab("Day") +ylab("Closing Price (US$)") +guides(colour=guide_legend(title="Forecast"))
aug %>%ACF(.resid) %>%autoplot() +ggtitle("ACF of residuals")
Se puede ver que los residuos no están autocorrelacionados (ACF), por lo tanto este método (NAIVE) es bueno para esta serie de tiempo, y los pronósticos derivados de este método pueden ser muy buenos.
google_2015 %>%model(MEAN(Close)) %>%gg_tsresiduals() +ggtitle("Diagnóstico de residuales para el modelo de Media")
Como podemos ver, el método de la media no es óptimo para esta serie de acciones, siendo que la media es 0, al graficar los residuales vemos que muestra un patón (igual al de los datos menos la media), existen claramente autocorrelaciones de una caminata aleatoria y no hay ninguna distribución normal. estos tres factores lo hacen un método nada efectivo para acciones.
Tests de Portmanteau de autocorrelación
Tests para analizar de una forma más formal la presencia o ausencia de autocorrelación en los residuales. Nos ayuda a ver si las primeras \(h\) autocorrelaciones son significativamente distintas de cero o no.
Test de Box-Pierce
Sumatoria de \(r^2_k\) donde \(h\) es el rezago máximo a considerar y \(T\) es la cantidad de observaciones en la muestra. Si \(r^2_k\) es pequeña, entonces el resultado será pequeño. Se recomienda usar \(h = 10\) para datos no estacionales, \(h = 2m\) para datos estaionales (m es el periodo estacional) o como máximo \(h = T/5\).
Hay que tomar en cuenta, que la prueba no es tan buena cuando \(h\) es grande (mayor a \(T/5\)).
Test de Ljung-Box
Generalmente más preciso que Box-Pierce, de igual forma, valores grandes en el resultado indica que no es ruido blanco y los residuales sí están autocorrelacionados.
Hipótesis nula
La hipótesis nula dice que la serie en cuestión no está autocorrelacionada, o sea, es riudo blanco. Para unos residuales, que la hipótesis nula sea cierta es bueno, ya que significa que no están correlacionados los residuales.
¿Cómo saber si es cierta o no?
Siendo \(\alpha\) el nivel de significancia o nivel máximo de error dispuestos a aceptar (usaremos 5%) y p-value un avlor resultante de las pruebas:
Si p-value es menor a \(\alpha\), entonces rechazamos la hiótesis nula, de lo contrario, la aceptamos.
O sea:
\(p-value > \alpha\), H0 es cierta, residuales no están autocorrelacionados.
\(p-value < \alpha\), H0 es falsa, residuales están autocorrelacionados.
Ejemplos:
# lag=h and fitdf=Kaug %>%features(.resid, box_pierce, lag=10, dof=0)
# A tibble: 1 × 4
Symbol .model bp_stat bp_pvalue
<chr> <chr> <dbl> <dbl>
1 GOOG NAIVE(Close) 7.74 0.654
aug %>%features(.resid, ljung_box, lag=10, dof=0)
# A tibble: 1 × 4
Symbol .model lb_stat lb_pvalue
<chr> <chr> <dbl> <dbl>
1 GOOG NAIVE(Close) 7.91 0.637
Como podemos ver para el método Naive, el p-value es mayor al 5%, por lo tanto sus residuos no están autocrrelacionados, son ruido blanco.
google_2015 %>%model(MEAN(Close)) %>%augment() %>%features(.resid, box_pierce, lag =10, dof =0)
# A tibble: 1 × 4
Symbol .model bp_stat bp_pvalue
<chr> <chr> <dbl> <dbl>
1 GOOG MEAN(Close) 2024. 0
google_2015 %>%model(MEAN(Close)) %>%augment() %>%features(.resid, ljung_box, lag =10, dof =0)
# A tibble: 1 × 4
Symbol .model lb_stat lb_pvalue
<chr> <chr> <dbl> <dbl>
1 GOOG MEAN(Close) 2083. 0
Para el caso de Google, no se logra distinguir los residuales del ruido blanco.
Intervalos de predicción
\(c\) es el porcentaje de cobertura de probabilidad. Normalmente se usará 80% y 95%.
Al tener incertidumbre en el ponóstico, entre maoyr sea el intervalo, menos preciso será nuestro pronóstico.
Intervalos de predicción de un paso
Cunado ser realizan pronósitocs de un paso, la desviación estándar del pronóstico es prácticamente la misma que la desviación estándar de los residuos.
Intervalos de predicción de paso múltiple
Entre más aumenta el horizonte de pronóstico, mayor será el intervalo de pronóstico, Esto es, \(\simga_h\) incrementa con \(h\), entonces necesitamos estimaciones de \(\sigma_h\).
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
of ggplot2 3.3.4.
Con esto, podemos obtener intervalos de predicción, al calcular los percentiles de los escenarios futuros. El resultado es el ntervalo de predicción bootstrapped. Esto se puede lograr con forecast:
fc %>%autoplot(google_2015) +ggtitle("Google closing stock price")
Pronósticos con transformaciones
Si se le hace un pronóstico a una serie anteriormente transformada (Box-cox por ejemplo), al terminar el pronóstico, se tiene que hacer una transformación inversa para ver la versión original del pronóstico.
Transformación inversa de Box-Cox:
\(y_t=\)
\({exp(w_t), \text{ si }\lambda = 0}.\)
\((\lambda w_t^\lambda+1)^{1/\lambda}, \text{ en otro caso}.\)
La paquetería fable realiza la transofrmación ivnersa en automático, cuando se especifica en la definición del modelo.
Ajustes por sesgo
Un problema al realizar transformaciones matemáticas, como Box–Cox es que las estimaciones puntuales re transformadas ya no presentan la media de la distribución de predicción, sino que la mediana. Esto puede llegar a ser problema dependiendo del contexto, pero normalmente no causa mucha diferencia.
Ajustar por sesgo es la transformación inversa de la media, entre más grande sea la varianza, mayor será la diferencia entre la media y la mediana.
eggs <-as_tsibble(fma::eggs)
Registered S3 method overwritten by 'quantmod':
method from
as.zoo.data.frame zoo
dcmp %>%model(NAIVE(season_adjust)) %>%forecast() %>%autoplot(dcmp) +ylab("New orders index") +ggtitle("Pronóstico naïve de los datos desestacionalizados")
Le podemos agregar nuevamente la estacionalidad con la función decomposition_model():
Usamos accuracy para ver los errores en el ajuste de modelos (mabble) a los datos de entrenamiento:
accuracy(beer_fit)
# A tibble: 4 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Mean Training 0 43.6 35.2 -0.937 7.89 2.46 2.60 -0.109
2 Naïve Training 4.76e- 1 65.3 54.7 -0.916 12.2 3.83 3.89 -0.241
3 Seasonal naïve Training -2.13e+ 0 16.8 14.3 -0.554 3.31 1 1 -0.288
4 Drift Training -5.41e-15 65.3 54.8 -1.03 12.2 3.83 3.89 -0.241
Ahora la usamos para ver los errores de pronóstico (fabble) comparado a los datos de prueba:
beer_fc |>accuracy(recent_production)
# A tibble: 4 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Drift Test -54.0 64.9 58.9 -13.6 14.6 4.12 3.87 -0.0741
2 Mean Test -13.8 38.4 34.8 -3.97 8.28 2.44 2.29 -0.0691
3 Naïve Test -51.4 62.7 57.4 -13.0 14.2 4.01 3.74 -0.0691
4 Seasonal naïve Test 5.2 14.3 13.4 1.15 3.17 0.937 0.853 0.132
Los errores de ajuste (datos de entrenamiento) sirven para ver que modelos utilizar para pronosticar.
Los errores de pronóstico (datos de prueba) nos da mayor claridad de cuáles modelos parecen producir los mejores pronósticos.
Después de esto se recalculan los modelos utilizando toda la información (sin separar en conjuntos de entrenamiento y prueba) y se producen los pronósticos reales a presentar.
